home *** CD-ROM | disk | FTP | other *** search
/ Chip 2007 January, February, March & April / Chip-Cover-CD-2007-02.iso / Pakiet bezpieczenstwa / mini Pentoo LiveCD 2006.1 / mpentoo-2006.1.iso / livecd.squashfs / usr / lib / python2.4 / site-packages / Numeric / FFT / FFT.py < prev    next >
Encoding:
Python Source  |  2006-01-20  |  10.4 KB  |  314 lines
  1. """
  2. Discrete Fourier Transforms - FFT.py 
  3.  
  4. The underlying code for these functions is an f2c translated and modified
  5. version of the FFTPACK routines.
  6.  
  7. fft(a, n=None, axis=-1) 
  8. inverse_fft(a, n=None, axis=-1) 
  9. real_fft(a, n=None, axis=-1) 
  10. inverse_real_fft(a, n=None, axis=-1)
  11. hermite_fft(a, n=None, axis=-1)
  12. inverse_hermite_fft(a, n=None, axis=-1)
  13. fftnd(a, s=None, axes=None)
  14. inverse_fftnd(a, s=None, axes=None)
  15. real_fftnd(a, s=None, axes=None)
  16. inverse_real_fftnd(a, s=None, axes=None)
  17. fft2d(a, s=None, axes=(-2,-1)) 
  18. inverse_fft2d(a, s=None, axes=(-2, -1))
  19. real_fft2d(a, s=None, axes=(-2,-1)) 
  20. inverse_real_fft2d(a, s=None, axes=(-2, -1))
  21. """
  22. import Numeric, fftpack, copy
  23.  
  24. _fft_cache = {}
  25. _real_fft_cache = {}
  26.  
  27. def _raw_fft(a, n=None, axis=-1, init_function=fftpack.cffti, 
  28.              work_function=fftpack.cfftf, fft_cache = _fft_cache ):
  29.     a = Numeric.asarray(a)
  30.  
  31.     if n == None: n = a.shape[axis]
  32.  
  33.     try:
  34.         wsave = fft_cache[n]
  35.     except(KeyError):
  36.         wsave = init_function(n)
  37.         fft_cache[n] = wsave
  38.  
  39.     if a.shape[axis] != n:
  40.         s = list(a.shape)
  41.         if s[axis] > n:
  42.             index = [slice(None)]*len(s)
  43.             index[axis] = slice(0,n)
  44.             a = a[index]
  45.         else:
  46.             index = [slice(None)]*len(s)
  47.             index[axis] = slice(0,s[axis])
  48.             s[axis] = n
  49.             z = Numeric.zeros(s, a.typecode())
  50.             z[index] = a
  51.             a = z
  52.  
  53.     if axis != -1:
  54.         a = Numeric.swapaxes(a, axis, -1)
  55.     r = work_function(a, wsave)
  56.     if axis != -1:
  57.         r = Numeric.swapaxes(r, axis, -1)
  58.     return r
  59.  
  60.  
  61. def fft(a, n=None, axis=-1):
  62.     """fft(a, n=None, axis=-1) 
  63.  
  64.     Will return the n point discrete Fourier transform of a. n defaults to the
  65.     length of a. If n is larger than a, then a will be zero-padded to make up
  66.     the difference. If n is smaller than a, the first n items in a will be
  67.     used.
  68.  
  69.     The packing of the result is "standard": If A = fft(a, n), then A[0]
  70.     contains the zero-frequency term, A[1:n/2+1] contains the
  71.     positive-frequency terms, and A[n/2+1:] contains the negative-frequency
  72.     terms, in order of decreasingly negative frequency. So for an 8-point
  73.     transform, the frequencies of the result are [ 0, 1, 2, 3, 4, -3, -2, -1].
  74.  
  75.     This is most efficient for n a power of two. This also stores a cache of
  76.     working memory for different sizes of fft's, so you could theoretically
  77.     run into memory problems if you call this too many times with too many
  78.     different n's."""
  79.  
  80.     return _raw_fft(a, n, axis, fftpack.cffti, fftpack.cfftf, _fft_cache)
  81.  
  82.  
  83. def inverse_fft(a, n=None, axis=-1):
  84.     """inverse_fft(a, n=None, axis=-1) 
  85.  
  86.     Will return the n point inverse discrete Fourier transform of a.  n
  87.     defaults to the length of a. If n is larger than a, then a will be
  88.     zero-padded to make up the difference. If n is smaller than a, then a will
  89.     be truncated to reduce its size.
  90.  
  91.     The input array is expected to be packed the same way as the output of
  92.     fft, as discussed in it's documentation.
  93.  
  94.     This is the inverse of fft: inverse_fft(fft(a)) == a within numerical
  95.     accuracy.
  96.  
  97.     This is most efficient for n a power of two. This also stores a cache of
  98.     working memory for different sizes of fft's, so you could theoretically
  99.     run into memory problems if you call this too many times with too many
  100.     different n's."""
  101.  
  102.     a = Numeric.asarray(a).astype(Numeric.Complex)
  103.     if n == None:
  104.         n = Numeric.shape(a)[axis]
  105.     return _raw_fft(a, n, axis, fftpack.cffti, fftpack.cfftb, _fft_cache) / n
  106.  
  107.  
  108. def real_fft(a, n=None, axis=-1):
  109.     """real_fft(a, n=None, axis=-1) 
  110.  
  111.     Will return the n point discrete Fourier transform of the real valued
  112.     array a. n defaults to the length of a. n is the length of the input, not
  113.     the output.
  114.  
  115.     The returned array will be the nonnegative frequency terms of the
  116.     Hermite-symmetric, complex transform of the real array. So for an 8-point
  117.     transform, the frequencies in the result are [ 0, 1, 2, 3, 4]. The first
  118.     term will be real, as will the last if n is even. The negative frequency
  119.     terms are not needed because they are the complex conjugates of the
  120.     positive frequency terms. (This is what I mean when I say
  121.     Hermite-symmetric.)
  122.  
  123.     This is most efficient for n a power of two."""
  124.  
  125.     a = Numeric.asarray(a).astype(Numeric.Float)
  126.     return _raw_fft(a, n, axis, fftpack.rffti, fftpack.rfftf, _real_fft_cache)
  127.  
  128.  
  129. def inverse_real_fft(a, n=None, axis=-1):
  130.     """inverse_real_fft(a, n=None, axis=-1)
  131.     
  132.     Will return the real valued n point inverse discrete Fourier transform of
  133.     a, where a contains the nonnegative frequency terms of a Hermite-symmetric
  134.     sequence. n is the length of the result, not the input. If n is not
  135.     supplied, the default is 2*(len(a)-1). If you want the length of the
  136.     result to be odd, you have to say so.
  137.  
  138.     If you specify an n such that a must be zero-padded or truncated, the
  139.     extra/removed values will be added/removed at high frequencies. One can
  140.     thus resample a series to m points via Fourier interpolation by: a_resamp
  141.     = inverse_real_fft(real_fft(a), m).
  142.  
  143.     This is the inverse of real_fft:
  144.     inverse_real_fft(real_fft(a), len(a)) == a
  145.     within numerical accuracy."""
  146.  
  147.     a = Numeric.asarray(a).astype(Numeric.Complex)
  148.     if n == None:
  149.         n = (Numeric.shape(a)[axis] - 1) * 2
  150.     return _raw_fft(a, n, axis, fftpack.rffti, fftpack.rfftb,
  151.                     _real_fft_cache) / n
  152.  
  153.  
  154. def hermite_fft(a, n=None, axis=-1):
  155.     """hermite_fft(a, n=None, axis=-1)
  156.     inverse_hermite_fft(a, n=None, axis=-1)
  157.  
  158.     These are a pair analogous to real_fft/inverse_real_fft, but for the
  159.     opposite case: here the signal is real in the frequency domain and has
  160.     Hermite symmetry in the time domain. So here it's hermite_fft for which
  161.     you must supply the length of the result if it is to be odd.
  162.  
  163.     inverse_hermite_fft(hermite_fft(a), len(a)) == a
  164.     within numerical accuracy."""
  165.  
  166.     a = Numeric.asarray(a).astype(Numeric.Complex)
  167.     if n == None:
  168.         n = (Numeric.shape(a)[axis] - 1) * 2
  169.     return inverse_real_fft(Numeric.conjugate(a), n, axis) * n
  170.  
  171.  
  172. def inverse_hermite_fft(a, n=None, axis=-1):
  173.     """hermite_fft(a, n=None, axis=-1)
  174.     inverse_hermite_fft(a, n=None, axis=-1)
  175.  
  176.     These are a pair analogous to real_fft/inverse_real_fft, but for the
  177.     opposite case: here the signal is real in the frequency domain and has
  178.     Hermite symmetry in the time domain. So here it's hermite_fft for which
  179.     you must supply the length of the result if it is to be odd.
  180.  
  181.     inverse_hermite_fft(hermite_fft(a), len(a)) == a
  182.     within numerical accuracy."""
  183.     
  184.     a = Numeric.asarray(a).astype(Numeric.Float)
  185.     if n == None:
  186.         n = Numeric.shape(a)[axis]
  187.     return Numeric.conjugate(real_fft(a, n, axis))/n
  188.  
  189.  
  190. def _cook_nd_args(a, s=None, axes=None, invreal=0):
  191.     if s is None:
  192.         shapeless = 1
  193.         if axes == None:
  194.             s = list(a.shape)
  195.         else:
  196.             s = Numeric.take(a.shape, axes)
  197.     else:
  198.         shapeless = 0
  199.     s = list(s)
  200.     if axes == None:
  201.         axes = range(-len(s), 0)
  202.     if len(s) != len(axes):
  203.         raise ValueError, "Shape and axes have different lengths."
  204.     if invreal and shapeless:
  205.         s[axes[-1]] = (s[axes[-1]] - 1) * 2
  206.     return s, axes
  207.  
  208.  
  209. def _raw_fftnd(a, s=None, axes=None, function=fft):
  210.     a = Numeric.asarray(a)
  211.     s, axes = _cook_nd_args(a, s, axes)
  212.     itl = range(len(axes))
  213.     itl.reverse()
  214.     for ii in itl:
  215.         a = function(a, n=s[ii], axis=axes[ii])
  216.     return a
  217.  
  218.  
  219. def fftnd(a, s=None, axes=None):
  220.     """fftnd(a, s=None, axes=None)
  221.  
  222.     The n-dimensional fft of a. s is a sequence giving the shape of the input
  223.     an result along the transformed axes, as n for fft. Results are packed
  224.     analogously to fft: the term for zero frequency in all axes is in the
  225.     low-order corner, while the term for the Nyquist frequency in all axes is
  226.     in the middle.
  227.  
  228.     If neither s nor axes is specified, the transform is taken along all
  229.     axes. If s is specified and axes is not, the last len(s) axes are used.
  230.     If axes are specified and s is not, the input shape along the specified
  231.     axes is used. If s and axes are both specified and are not the same
  232.     length, an exception is raised."""
  233.  
  234.     return _raw_fftnd(a,s,axes,fft)
  235.  
  236.  
  237. def inverse_fftnd(a, s=None, axes=None):
  238.     """inverse_fftnd(a, s=None, axes=None)
  239.     
  240.     The inverse of fftnd."""
  241.     
  242.     return _raw_fftnd(a, s, axes, inverse_fft)
  243.  
  244.  
  245. def fft2d(a, s=None, axes=(-2,-1)):
  246.     """fft2d(a, s=None, axes=(-2,-1)) 
  247.     
  248.     The 2d fft of a. This is really just fftnd with different default
  249.     behavior."""
  250.  
  251.     return _raw_fftnd(a,s,axes,fft)
  252.  
  253.  
  254. def inverse_fft2d(a, s=None, axes=(-2,-1)):
  255.     """inverse_fft2d(a, s=None, axes=(-2, -1))
  256.  
  257.     The inverse of fft2d. This is really just inverse_fftnd with different
  258.     default behavior."""
  259.  
  260.     return _raw_fftnd(a, s, axes, inverse_fft)
  261.  
  262.  
  263. def real_fftnd(a, s=None, axes=None):
  264.     """real_fftnd(a, s=None, axes=None)
  265.  
  266.     The n-dimensional discrete Fourier transform of a real array a. A real
  267.     transform as real_fft is performed along the axis specified by the last
  268.     element of axes, then complex transforms as fft are performed along the
  269.     other axes."""
  270.     
  271.     a = Numeric.asarray(a).astype(Numeric.Float)
  272.     s, axes = _cook_nd_args(a, s, axes)
  273.     a = real_fft(a, s[-1], axes[-1])
  274.     for ii in range(len(axes)-1):
  275.         a = fft(a, s[ii], axes[ii])
  276.     return a
  277.  
  278.  
  279. def real_fft2d(a, s=None, axes=(-2,-1)):
  280.     """real_fft2d(a, s=None, axes=(-2,-1)) 
  281.  
  282.     The 2d fft of the real valued array a. This is really just real_fftnd with
  283.     different default behavior."""
  284.     
  285.     return real_fftnd(a, s, axes)
  286.  
  287.  
  288. def inverse_real_fftnd(a, s=None, axes=None):
  289.     """inverse_real_fftnd(a, s=None, axes=None)
  290.  
  291.     The inverse of real_fftnd. The transform implemented in inverse_fft is
  292.     applied along all axes but the last, then the transform implemented in
  293.     inverse_real_fft is performed along the last axis. As with
  294.     inverse_real_fft, the length of the result along that axis must be
  295.     specified if it is to be odd."""
  296.     
  297.     a = Numeric.asarray(a).astype(Numeric.Complex)
  298.     s, axes = _cook_nd_args(a, s, axes, invreal=1)
  299.     for ii in range(len(axes)-1):
  300.         a = inverse_fft(a, s[ii], axes[ii])
  301.     a = inverse_real_fft(a, s[-1], axes[-1])
  302.     return a
  303.  
  304.  
  305. def inverse_real_fft2d(a, s=None, axes=(-2,-1)):
  306.     """inverse_real_fft2d(a, s=None, axes=(-2, -1))
  307.  
  308.     The inverse of real_fft2d. This is really just inverse_real_fftnd with
  309.     different default behavior."""
  310.     
  311.     return inverse_real_fftnd(a, s, axes)
  312.  
  313.  
  314.